{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Chapter 16\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Code from previous notebooks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_func(state, t, system):\n",
    "    \"\"\"Update the thermal transfer model.\n",
    "    \n",
    "    state: State (temp)\n",
    "    t: time\n",
    "    system: System object\n",
    "    \n",
    "    returns: State (temp)\n",
    "    \"\"\"\n",
    "    r, T_env, dt = system.r, system.T_env, system.dt\n",
    "    \n",
    "    T = state.T\n",
    "    T += -r * (T - T_env) * dt\n",
    "    \n",
    "    return State(T=T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_simulation(system, update_func):\n",
    "    \"\"\"Runs a simulation of the system.\n",
    "    \n",
    "    Add a TimeFrame to the System: results\n",
    "    \n",
    "    system: System object\n",
    "    update_func: function that updates state\n",
    "    \"\"\"\n",
    "    init = system.init\n",
    "    t_0, t_end, dt = system.t_0, system.t_end, system.dt\n",
    "    \n",
    "    frame = TimeFrame(columns=init.index)\n",
    "    frame.row[t_0] = init\n",
    "    ts = linrange(t_0, t_end, dt)\n",
    "    \n",
    "    for t in ts:\n",
    "        frame.row[t+dt] = update_func(frame.row[t], t, system)\n",
    "    \n",
    "    return frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_system(T_init, r, volume, t_end):\n",
    "    \"\"\"Makes a System object with the given parameters.\n",
    "\n",
    "    T_init: initial temperature in degC\n",
    "    r: heat transfer rate, in 1/min\n",
    "    volume: volume of liquid in mL\n",
    "    t_end: end time of simulation\n",
    "    \n",
    "    returns: System object\n",
    "    \"\"\"\n",
    "    init = State(T=T_init)\n",
    "                   \n",
    "    return System(init=init,\n",
    "                  r=r, \n",
    "                  volume=volume,\n",
    "                  temp=T_init,\n",
    "                  t_0=0, \n",
    "                  t_end=t_end, \n",
    "                  dt=1,\n",
    "                  T_env=22)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using `root_bisect`\n",
    "\n",
    "As a simple example, let's find the roots of this function; that is, the values of `x` that make the result 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func(x):\n",
    "    return (x-1) * (x-2) * (x-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`modsim.py` provides `root_bisect`, which searches for a root by bisection.  The first argument is the function whose roots we want.  The second argument is an interval that contains a root."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>converged</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>root</th>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "converged    True\n",
       "root            1\n",
       "dtype: object"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res = root_bisect(func, [0.5, 1.5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is an object that contains the root that was found and other information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.root"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we provide a different interval, we find a different root."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>converged</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>root</th>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "converged    True\n",
       "root            2\n",
       "dtype: object"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res = root_bisect(func, [1.5, 2.5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.root"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the interval doesn't contain a root, the results explain the error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>converged</th>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>flag</th>\n",
       "      <td>4.000000 and 5.000000 do not bracket a root</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "converged                                          False\n",
       "flag         4.000000 and 5.000000 do not bracket a root\n",
       "dtype: object"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res = root_bisect(func, [4, 5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to find the value of `r` that makes the final temperature 70, so we define an \"error function\" that takes `r` as a parameter and returns the difference between the final temperature and the goal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def error_func1(r):\n",
    "    \"\"\"Runs a simulation and returns the `error`.\n",
    "    \n",
    "    r: heat transfer rate, in 1/min\n",
    "    \n",
    "    returns: difference between final temp and 70 C\n",
    "    \"\"\"\n",
    "    system = make_system(T_init=90, r=r, volume=300, t_end=30)\n",
    "    results = run_simulation(system, update_func)\n",
    "    T_final = get_last_value(results.T)\n",
    "    return T_final - 70"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With `r=0.01`, we end up a little too warm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.2996253904030937"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "error_func1(r=0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With `r=0.02`, we end up too cold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-10.907066281994297"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "error_func1(r=0.02)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The return value from `root_bisect` is an array with a single element, the estimated value of `r`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>converged</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>root</th>\n",
       "      <td>0.0115431</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "converged         True\n",
       "root         0.0115431\n",
       "dtype: object"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res = root_bisect(error_func1, [0.01, 0.02])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.011543084681034089"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r_coffee = res.root"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we run the simulation with the estimated value of `r`, the final temperature is 70 C, as expected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "69.99999985860761"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coffee = make_system(T_init=90, r=r_coffee, volume=300, t_end=30)\n",
    "results = run_simulation(coffee, update_func)\n",
    "T_final = get_last_value(results.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  When you call `root_bisect`, it calls `error_func1` several times.  To see how this works, add a print statement to `error_func1` and run `root_bisect` again."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Repeat this process to estimate `r_milk`, given that it starts at 5 C and reaches 20 C after 15 minutes.  \n",
    "\n",
    "Before you use `root_bisect`, you might want to try a few values for `r_milk` and see how close you can get by trial and error.  Here's an initial guess to get you started:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "18.499850754390966"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r_milk = 0.1\n",
    "milk = make_system(T_init=5, r=r_milk, volume=50, t_end=15)\n",
    "results = run_simulation(milk, update_func)\n",
    "T_final = get_last_value(results.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "def error_func2(r):\n",
    "    \"\"\"Runs a simulation and returns the `error`.\n",
    "    \n",
    "    r: heat transfer rate, in 1/min\n",
    "    \n",
    "    returns: difference between final temp and 20C\n",
    "    \"\"\"\n",
    "    system = make_system(T_init=5, r=r, volume=50, t_end=15)\n",
    "    results = run_simulation(system, update_func)\n",
    "    T_final = get_last_value(results.T)\n",
    "    return T_final - 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.500149245609034"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "error_func2(r=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.4018656744898585"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "error_func2(r=0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>converged</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>root</th>\n",
       "      <td>0.132961</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "converged        True\n",
       "root         0.132961\n",
       "dtype: object"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "res = root_bisect(error_func2, [0.1, 0.2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.13296079039573666"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "r_milk = res.root"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "20.00000003602163"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "milk = make_system(T_init=5, r=r_milk, volume=50, t_end=15)\n",
    "results = run_simulation(milk, update_func)\n",
    "T_final = get_last_value(results.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mixing liquids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function takes `System` objects that represent two liquids, computes the temperature of the mixture, and returns a new `System` object that represents the mixture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mix(system1, system2):\n",
    "    \"\"\"Simulates the mixture of two liquids.\n",
    "    \n",
    "    system1: System representing coffee\n",
    "    system2: System representing milk\n",
    "    \n",
    "    returns: System representing the mixture\n",
    "    \"\"\"\n",
    "    assert system1.t_end == system2.t_end\n",
    "    \n",
    "    V1, V2 = system1.volume, system2.volume\n",
    "    T1, T2 = system1.temp, system2.temp\n",
    "    \n",
    "    V_mix = V1 + V2\n",
    "    T_mix = (V1 * T1 + V2 * T2) / V_mix\n",
    "    \n",
    "    return make_system(T_init=T_mix,\n",
    "                       r=system1.r,\n",
    "                       volume=V_mix,\n",
    "                       t_end=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`mix` requires the `System` objects to have `temp` as a system variable.  `make_system` initializes this variable;\n",
    "the following function makes sure it gets updated when we run a simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_and_set(system):\n",
    "    \"\"\"Run a simulation and set the final temperature.\n",
    "    \n",
    "    system: System\n",
    "    \n",
    "    returns: TimeFrame\n",
    "    \"\"\"\n",
    "    results = run_simulation(system, update_func)\n",
    "    system.temp = get_last_value(results.T)\n",
    "    return results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mixing immediately\n",
    "\n",
    "Next here's what we get if we add the milk immediately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "90"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coffee = make_system(T_init=90, r=r_coffee, volume=300, t_end=30)\n",
    "coffee.temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "milk = make_system(T_init=5, r=r_milk, volume=50, t_end=30)\n",
    "milk.temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "77.85714285714286"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mix_first = mix(coffee, milk)\n",
    "mix_first.temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "61.4285713124277"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mix_results = run_and_set(mix_first)\n",
    "mix_first.temp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mixing at the end\n",
    "\n",
    "First we'll see what happens if we add the milk at the end.  We'll simulate the coffee and the milk separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "69.99999985860761"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coffee_results = run_and_set(coffee)\n",
    "coffee.temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "21.76470589082862"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "milk_results = run_and_set(milk)\n",
    "milk.temp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the results look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving figure to file figs/chap16-fig01.pdf\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZhcVYH38W8tXb3ve/YActghYQkIEWQRERdU3kGJgnsQVERxGRiDyCjIpqIoRHnZlGFwgIFHXxFhQIEBjBBkCTkBQkjS6X3fu2t5/zi3qqs73Z1K0l29/T7PU09V3Xvr1rl9Sf045557ji8WiyEiIjLd+Ke6ACIiIqNRQImIyLSkgBIRkWlJASUiItNScKoLsLeMMZnA0UAtEJni4oiIyO4JANXAOmttf/KKGR9QuHB6aqoLISIie2Ul8HTygtkQULUAv/vd76iqqprqsoiIyG6oq6tj1apV4P2WJ5sNARUBqKqqYsGCBVNdFhER2TM7XaJRJwkREZmWFFAiIjItKaBERGRaUkCJiMi0pIASEZFpSQElIiLTkgIK6Ood5PW3W2jv6t/1xiIikhaz4T6ovfbT/3iR51+rA6AoL5NFVfksri5gcVU+iyoLWFSVT252xhSXUkRkblFAAScftZDWzj621XfS1tVP25v9vPxm07BtyoqyXWBVecFVlc/CinyyMvUnFBGZDPp1Bd592Dzefdg8otEYTW29vFPXwTt1nbxT18HWuk621XfS1NZLU1svL2xsSHzO54PKkpxELWtRVT6LKvNZUJlPZkZgCo9IRGTmU0Al8ft9VJTkUFGSw9EHDY3rF4nGqGvu5p3aDrbWdyaeaxq6qGvuoa65h79vqBvajw+qSnO90CpgUaULr/nleYQUXCIiKVFApSDg9zG/PI/55Xm8O2n5YDjKjqYuttZ1uke9q3HtaOpOPJ57dXhwVZflsqiqgIWV+QouEZFxKKD2QkbQz+KqAhZXFQxbPhiOsL3BC676TrbWdbCtvpPapm5qGt3j2VeGBu5VcImI7EwBNQkyggGWzitk6bzCYcsHBiPUNO5ecFWW5iYCa2GleyyoyCMrpFMnIrObfuXSKJQxfnC9U9fJ9vp4eHVS29xNbZN7xLvBw1DnjIWVrifhQi/AFlTkkZOl7vAiMjsooKaBsYJrMBxhR2N3IrC2eeG1o3Goc8a6DfXDPlNWlM3CijwWej0K47Wu/JxQOg9JRGSvKaCmsYxgwN0wXD38Glc4EmVHYxfbGrrYVt/JNq/JsKaxK9Edfv2mxmGfKcrPZGFFPgsq81hYEe8On0dJQRY+ny+dhyUikhIF1AwUDPhd9/URnTMi0Rj1Ld2JwNre0OWe6ztp6+ynrbOfV94afgNyTlYw0Uy4sDKPBV6zYUVJDgG/gktEpo4CahYJ+H3MK8tjXlkeKw6pTiyP34C8raGTbfVdbG9wzYXb6jvp7BnEbm3Fbm0dtq9Q0M+88jwWVOQlrnUtqMxjXnmebkIWkbRQQM0ByTcgH3lAZWJ5LBajvWuAbQ1DnTO213exraGT5vY+ttR2sKW2Y9i+4h00FlTkJ8Ir/qzrXCIykRRQc5jP56MoP5Oi/EwO3bds2LqevkG2N8RrW+5a1/aGTmq9zhl1zT384/XhHTQK80KJ4FpQ4TUZVuRTXpSNX82FIrKbFFAyqpysDPZfVMz+i4qHLR8MR6ltch00tte78KppdNe72rsGaO9q5rXNzcM+E8oIsMBrLoyH1/yKPOaV5+p+LhEZk34dZLdkBEfvoBGNxmhu73M1rgYXWDVeL8PWzn4272hn8472nfZXUZydqHXNTwqw4vxM9S4UmeMUUDIh/H4f5cXZlBdns8xUDFvX1TtITUO8tuWaDbc3dFHb1E1Day8Nrb28aBuGfSYnK+jGP4yHVrkLseqyXA3/JDJHKKBk0uVlZ2AWl2AWlwxbHo5EqW/pYbvXJT5+zWt7QxddvYO8sa2NN7a1DfuMzwcVxTlDNa7yoSZD1bpEZhcFlEyZYMCfGCV+RdLyeO9CV9saqnXVNHRR19JDvfdInpsLIDszyPyKPOaXDYWXrnWJzFz6VyvTTnLvwoP3KR22bjAcpa65O1Hbqml017pqGrvo7BnkzW1tvDmi1gVuCKh4YMVDcX5FHmVF2bohWWSaUkDJjJIR9CfGF4TqYevau/qHBVa89lXb1J0YAuqlNxp32t+8stzETcnJ4aX7ukSmVloDyhhzLHATYIBG4Bpr7W+MMSHgF8DZQAS40Vp7dTrLJjNfYV4mhXmZHLR0eK0rEr/WlRReNY1d7GjsoqWjn3fqOnmnrnOn/eXnhFjgNRHOL3ejaMwvdx01NJqGyORLW0AZY/zAQ8A3rbW/NcYcDTxljFkHfAIXWvsChcAjxpgaa+1d6SqfzF6BgBu2aV55Hhw0fF1P3yA7GrvZ7gVWTUMXNU3udWfPAK9vaeH1LS077bO8OJv5ZTuHV0VxNoGAP01HJjK7pbMGVQxUAD5jjA+IAWFgADgf+Iy1thVoNcZcD6wGFFAyqXKyMthvYRH7LSwatjwWi9HS0ceOxu6kGpd7XdfcTWNrL42tOzcZBgM+KktyhwXXvDL3uqQgSyNqiOyGtAWUtbbZGPML4E7gdiAAXAzU4i4mbEjafCNwaLrKJjKSz+ejtDCb0sJsDt1v+DBQkUiU+taeEeHlal9N7X2JZesYPhRUKCPgXe/ywqssl2qvFlaUpy7yIiOlu4mvDzgXuB94N/AAEO9y1ZO0eQ+Qk66yieyOQMCfGDX+qAMrh63rGwhT19yTCK14iO1ockNBjTYAL7gu8vPKc7395npNku59Qa46a8jclM4mvo8Bx1trv+W9/6sx5jZc8x5AdtK2OUBXGssmMiGyQkGWVBewZMQkk+BG1NgRD66mbnY0drPDu97V3Rfmre3tvLV95+Gg8rIzmFeeS3VpPLRyE02HeeppKLNYOgNqIZA5YlkY15uvDtdJosZbfgDDm/xEZry87NEH4I3FYnR0DwwFVlN3IsRqm9yoGpu2trFp6873d+XnZDCvLI9qr7ZVXeYFmMJLZoF0BtSjwNXGmC8BvwaWA18EvgBsBa4wxrwM5AGXAj9LY9lEpozP50t0kT9w6fDhoGKxGG2d/SNCywVZbVP3mBNOgusm765z5Saeq73rXmo2lJkgnZ0kXjPGfAy4CrgOV2v6rrX2IWPMn4EbgNcAP7AWuCVdZROZrnw+H8UFWRQXZO00qka8p6ELrOEBVtvcTWfPAHbrwKjhlZedkRRYXoCVuhpYYV5IHTZkWvDFYrGpLsNeMcYsAd5+/PHHWbBgwVQXR2RaiMVitHb2u/Bq7KK22YVYbWM3tc1d9PZHxvxsdmaQ6tLcYQFWXZZLdWmuusrLhNu+fTunnHIKwFJr7ZbkdRrqSGQW8vl8lBRkUTJGzautq58djd3UNXs1rqZudjR3U+t12Bhr/q5Q0E9lqatxVZUOBVd1Wa5uUpYJp4ASmWN8Ph/F+VkU548eXp09g9TFa1xeR40dTS7M2rsG2Fbfybb6nYeG8vt9VBbnUFWak6h1VZW6AKsszdGI8rLb9F+MiCT4fD4KckMU5IZ26m0Ibmio2qZu6pp72NHURV1zT+KaV1NbL7XN7vX6TY07fbakIJOq0qGalwuvHKpKcynI1XUv2ZkCSkRSlpOVwb4Lith3QdFO6wYGI9S3DAVWXfy5uZv6lh5aOvpp6ehnw9s7j22YkxWkqiSXqrIcr8Y1FF7lRWo6nKsUUCIyIUIZgaSpUIaLRGM0ezWs+HWvuuYeV+Nq6qZnnOteAb+PiuIcKktdeFV5wVXlvc7JykjH4ckUUECJyKQL+H1UlORQUZLD4e8qH7YufqNyXbMLrbrmeM3LvW5u70s0Hb7Ezk2H+TkhqstyqCpx17qSr3uVFmpCyplMASUiUyr5RmWzuGSn9f2DERpaehLNhnUtPYkwq/fu9+rcOjDqSBvBgKt9VXmBVVWSQ2VpLlUlbllutmpf05kCSkSmtcxxmg6j0RitnX0urFqGmg3rvefEKBxN3aPuOz8nIxFYlV5oxZ/Li7MJ6trXlFJAiciM5fcPTYsysss8QF9/mPrWHuq95sK6Fu+1F2adPYN09rTx5rada19+H5QWZbumwxLXfT45xIryNUXKZFNAicislZUZZHFVAYurdh5dPj7OYZ0XWPVe06F77qG5vTcxMeUrb+2871BGgMqSbCq9ABv2KM0lT82He00BJSJzUvI4hyMH6QUYDEdobO11ta4Wd72rznuub3G1r231XWyrH31moNzsjGGhFb/+VVGcTUWJblxOhf5CIiKjyAgGvIkj80Zd39M3mKhtDQsw79HdO8jmmnY21+zcdR6gKD/ThZfXhb6yJCfRnb68KIeMoK5/KaBERPZATlYGS+cVsnRe4U7r4uMd1nvXvBpaexKv61t7aGztoa2zn7bOfuw7O4827/NBaUEWFV7tqyIpyCqKcygrmhsdOBRQIiITLHm8wwNG6ToficZoae+jvqXbhVdzT6L21djaQ1NbL03tfTS194068obf76Os0AVYRXG89pWdeD9bAkwBJSKSZgG/j/LibMqLs0ddH45EaWrrpb6lh4YWV+uKv25o6aG5o4+G1l4aWnuB5p0+H++BOBReOwfYTGhCVECJiEwzwYA/MZzTaOIdOFzToXtuaB0eYPEeiK9t3jnA4k2I5V6AlRdne8/e+6JsQhmByT7MXVJAiYjMMLvqwDEYdjWwYbWv1h4aW12trKV9qAnx9S07NyGC68RRUZzt1b5cDaw8qTaWjjEQFVAiIrNMRtCfmJNrNPEmxHhgNbbGO2+4901tvYlOHKMNIQWQl51BRXEOHzxhKaetWDwpx6GAEhGZY5KbEA8dZX28E0dD6/DwamjpoaG1l8bWHrp6B+nqbedv62sUUCIikh7DO3HsPIRUvBt9c1sf88pHr6VNBAWUiIjsluRu9JNp+vczFBGROSmlGpQx5gzgDOAooAKIAHXAOuAP1tonJ6uAIiIyN40bUMaYTwLfB0qAx4A/4e4KCwBlwOHAfcaYRuCH1tp7JrW0IiIyZ4wZUMaYPwPdwAXAX6210TG28+NqV182xnzWWnvapJRURETmlPFqUGustc/vagdecP0R+KMx5tgJK5mIiMxpY3aSSCWcRvnMc3tXHBEREWdX16D2Ba4ELrfWvpO0/A4gC/hO8nIREZGJMmYNyhjzLuBZXEeIkXdiPQEcCDxvjFk6ecUTEZG5arz7oK7CBdQya+2G5BXW2juBo4FXgB9MXvFERGSuGi+g3gP8wFobHm2ltXYAF07vnYyCiYjI3DZeQOUDO89FPNx2YOf5jkVERPbSeAH1OrCrbuPHAVsmrDQiIiKe8QLqN8CPjDFLRltpjNkHuAa4axLKJSIic9yY3cyttWuNMScBG4wx/4Ubd68dKMZ1kPgYbuijG9NQThERmWPGvQ/KWnuuMeZc4HPAFUARbiy+54FPWWsfmPwiiojIXLTL0cy9AWA1CKyIiKTVeDfqXm6MyU51R8aYXGPMmokploiIzHXjdZLoB141xvzIGHPUWBsZY5YZY27A9frrm+gCiojI3DReJ4nrjTG/B74DPGmM6QU2AE2ADygHDgFCuJ5877HWbhnvy4wx1cCvcDf39gFrrbXfM8aEgF8AZ+MmQ7zRWnv1Xh6biIjMYLvqJPEOcKEx5jvAycCRQCUQBf4J/Bh43Frbm+L3PQS84O2jGvirMeZ14FDAAPvibvx9xBhTY61VF3YRkTkqpSnfrbWduHB5aE+/yBizAtgHON5aOwi87XVj7wWuBz5jrW0FWo0x1wOr0T1WIiJz1njXoCbakbjBZb9vjKkxxrwFfBQXUNW45sO4jbhalYiIzFEp1aAmSAmwEvgrriZ1APAI0Oit70natgfISWPZRERkmklnQPUDHdba73vv/2mM+Q1wvvc+uUt7DtCVxrKJiMg0k84mvo1AjtdjLy6IGzG9DtdJIu4Ahjf5iYjIHJNyDcq7F+obwP7AWcA5wGZr7YMp7uIvuOa8G4wx38QF0ueBLwObgSuMMS8DecClwM9SLZuIiMw+KdWgjDGn464d9QIH4+59ygH+0xjzmVT2Ya3tA07EXX+qxV1/utZaez+wBngVeA03KO39wC27cyAiIjK7pFqD+nfgm9baW4wx/wJgrb3KGNOIu5H3jlR2Yq3dDJw5yvI+4CLvISIikvI1qIOAR0dZ/iiwZMJKIyIi4kk1oGqAZaMsPxnYOnHFERERcVJt4rsaWOvNohsA3u/NtHsRcMkklU1EROawlGpQ1trbgfOAj+AGeb0GOAk431q7dtJKJyIic1ZKNShjzLeB31lrT5jk8oiIiACpX4O6HMiazIKIiIgkSzWg/gh81RhTPJmFERERiUu1k8R+wCdwIdWBu2E3wVo7b6ILJiIic1uqAXWr9xAREUmLVCcsvG2yCyJ77uqrr+a+++6joqKCu+++m9WrV7NlyxY++9nP8rWvfW2qiyciskdS7cV3z3jrrbXnTkxxZE/89re/5bbbbuPYY4/l4YcfprOzk3Xr1hEMpnM2FRGRiZVqJ4nIiIcP2Bc3qrlGkphA69ev55xzzmHZsmWcfvrpPProo4TDYW666SZOPPFEVqxYwQUXXMD27duJRCIsW7aMcDjM6tWrOfnkk7n88supqanh6KOP5u2336auro6LLrqIFStWcOqpp3LHHXckvisSiXDLLbdwyimnsGLFCi6++GJaWlqm7uBFRJKk2sT36dGWG2MuwwXVjHHlb57jH6/Xp+37jjqwkiu+cGxK27a0tPDFL36RSy65hHPOOYd169axevVqzjzzTF555RV++9vfUlFRwXXXXccFF1zAgw8+yPr16zHGcO+993LggQfywAMPcOedd/LQQw8RiUT4+Mc/zlFHHcXf/vY3amtrWb16NUVFRZx11lncddddPPzww9x+++1UVFTw4x//mEsuuYQ777xzkv8qIiK7trcTFt4D/J+JKIjAE088QWVlJatWrSIYDHLcccdxzz338Nhjj3HhhReycOFCMjMz+fa3v82OHTt4+eWXx93fq6++ytatW/nud79LZmYmS5Ys4bOf/Sz33nsvAPfddx9f+cpXWLRoEVlZWXzrW99i3bp1bNmyJQ1HKyIyvr29SPEZoHMCypE2qdZmpkJzczPV1dXDlh1yyCH09vYyb95QT/5QKERFRQV1dXXj7q+mpobe3l6OPXbomKPRKEVFRQDs2LGDyy+/nDVr1iTWB4NBampqWLJkyQQckYjInku1k0QtEBuxOM97fHOiCzVXVVRUUF8/vPnxjjvuIBaLUVNTwxFHHAHAwMAA9fX1lJaW7nJ/paWlPP3004llLS0t9PX1JdavWbOGlStXJtZv2rRJ4SQi00KqTXz/Bnwv6fFvwIXAgdban0xS2eacE088kfr6en7/+98TiUR49tlnuemmm7jwwgv51a9+xbZt2+jv7+faa6+luLiY5cuXj7u/ww47jLy8PH75y18yMDBAS0sLF154ITfddBMAZ511FjfffDO1tbVEIhHWrl3LqlWrEgEmIjKVUm3iqwR+Yq0dNoKEMabAGHOttfbbE1+0uae4uJi1a9dy9dVXc80111BZWckNN9zACSecwODgIOeddx7t7e0sX76c22+/nVAoNO7+QqEQa9eu5Uc/+hErV67E5/Nx6qmnctlllwHwpS99iXA4zKpVq2hra2P//ffntttuo6CgIB2HKyIyLl8sNrLlzjHGVAP53tvXgROA5hGbLQPusNZmT1oJd8Gbl+rtxx9/nAULFkxVMUREZA9s376dU045BWCptXZL8rrxalDvBn7P0LWnZ5LWxXD3QgHcMSGlFBERSTJmQFlr7zfG7Ie7TrUJOA5oStokBnRZaxsmt4giIjIXjXsNylq7GcAYk2GtjYy2jbducDIKJyIic1eqnSRKjDHfBQ4CAt4yH5AJHAKUTELZRERkDku1m/la4F+At4H3Am8AUVzHiR9OTtFERGQuSzWgTgbOt9ZeCLwG/F9r7enAtbjOFCIiIhMq1YDKxNWawHU5P9J7/X+B4ye6UCIiIqkG1CaGguh1YIX3Osd7iIiITKhUO0lcD9xhjAkA/wn80xgTA44FnpqswomIyNyVUg3KWnsXcCrwurXWAh/B9dx7Fvj85BVPRrNjxw6WLVtGZ2cn27dvxxhDR0fHsNciIjNdqqOZPwR8x1q7EcBa+2fgz5NZMBnbvHnzWL9+PQDt7e1TXBoRkcmR6jWoE4CBySyIuDGpjjrqKO6++26OP/54jjnmGO6++25+97vfsXLlSlasWMEdd9yRUk0pEolwySWXsGrVKrq7u9N4FCIiEyPVa1A3AncaY24ENgPDRjW31m6a6IJNps0//PiY68rOWE3B8vcB0PHiozT96dYxt93n8vsTr7ff9i0G6jaPu00qOjs7efXVV3niiSd47LHHuPTSS/nIRz7C448/zlNPPcVXv/pVjjnmmHH3EY1Gueyyy2hsbOTXv/41OTnqxyIiM0+qAXWV9xzvyRcfQNbnvQ7s9AnZYxdeeCGhUIjjjjuOSCTCeeedRygU4r3vfS+RSAS/f/yK75VXXsmLL77In/70J4WTiMxYqQbUuya1FGmWaq2mYPn7ErWpXVnw+ev2pkjDxKdkDwRc7ufnu1lP4sE01hQpcbW1tXR2drJu3TpOPPHECSuXiEg6pdqL7y1r7Vu4mtJBQC0wmLRcJpDP59v1RuO4+eab+frXv86aNWvo6uqaoFKJiKRXSgFljMn3evJtBB4EqoCfGWPWG2PmTWYBZfdlZGTwqU99ivLycq699tqpLo6IyB5JtRff9UAhsIShDhLfALqAn058sWRv+f1+rrrqKh544AGee+65qS6OiMhuG3PK92TGmBrgw9baF4wxncDh1trNxphlwGPW2tLJLug4ZVuCpnwXEZmR9nTK92T5QM8Y61LdBwDGmCLgZWCNtfYOY0wI+AVwNhABbrTWXr07+xQRkdkn1Sa+PwNrvLH4AGLGmFLgOuAvu/mdtwDzk95fCRhgX+Bo4HxjzHm7uU8REZllUg2orwFLgSbc6OWPAFuBMuDrqX6ZMeZ8oAB4JWnx+cAPrbWtXvXuemB1qvsUEZHZKaXmOWttLXCsMeY0XDfzIG7ajUestdFU9mGMWQpcgZvg8BFvWRFQDWxI2nQjcGiqByAiIrNTqjWouM3AFsACr+5GOAWA3wKXWmvrklblec/J17d60BxTIiJzXqqjmRcB9wGn4Doy+AC/MeZB4Dxr7VgdKOK+B1hr7QMjlsdHMc1OWpaD674uIiJzWKo1qFtx15uOAbJwU8Afh+vYkMp9UJ8AzjbGtBlj2nBNeL8EfgjU4TpJxB3A8CY/ERGZg1LtIn4GcJK19sWkZX83xqzGXU/60ngfttYekPzeGPMS8FOvm3kXcIUx5mVck9+lwM9SPQAREZmdUq1BNQHFoywPMtRMt6fWAK8CrwHrgPtxXdFFRGQOS7UG9a/AWmPMvwJPA2FgOfAT4OfGmP3jG6YyN5S19oik133ARd5DREQESD2g/sN7vpfhc0EBXANcjeaGEhGRCTQn54MSEZHpL9UbdTXnk4iIpFWq90EdBtwIHIzrYj6MtbZkgsslIiJzXKpNfHfieuv9G9A3ecURERFxUg2o/YGjrbW6gVZERNIi1fugnsE174mIiKRFqjWoLwLPGGPOwA0YO2yQWGvtjya6YCIiMrelGlDfw02LcTywbMS6GKCAEhGRCZVqQH0C+Ji19qHJLIyIiEhcqtegWoFdDmEkIiIyUVKtQX0HuMkY803gLWAweaW1dmCiCyYiInNbqgF1A1AKrB9jvcbfExGRCZVqQH1qUkshIiIyQqpj8T0ef22MyQe6rLWxcT4iIiKyV1KtQWGM+QbwLaAc2N8Y8z2gHbjUWhuepPKJiMgclVIvPq9zxMXAt4F+b/FDwDnoHigREZkEqXYz/yJwgbX2brxRJKy1/w2cD5w7SWUTEZE5LNWAWgRsHGX5FkBTbYiIyIRLNaDWA/+S9D7eQeICxu56LiIissdS7SRxKfCIMeY9uAkLrzTGHAAcBJwxWYUTEZG5K6UalLX2WdycUOuBPwLFwJPAgdbapyatdCIiMmeNWYMyxqwBrrfW9gBYa+txM+qKiIhMuvFqUFcAeekqiIiISLLxAsqXtlKIiIiMsKtOEguMMVm72om1dusElUdERATYdUCt28V6H67LuUYzFxGRCbWrgHov0JyOgoiIiCQbL6BiwEZrbUO6CiMiIhKnThIiIjItjRdQdwK96SqIiIhIsjGb+Ky1n01nQURERJKlOlisiIhIWimgRERkWlJAiYjItKSAEhGRaUkBJSIi01KqExaKiMgYYrEYxKLEohGIRohFIhCLEcjJT2wz2NZALDzgbReFqLd9LEogv5iMwgoAIt3t9Ne+6baJRYnF3Lbxz+WaFfhDbojUnrfWM9haB7EoxGJDn4lGySiqIO/gEwCIDvTR9sz9QHwbr7zec/7hp5BZtRSA7jf+QffG54H4djGvnDH8oSzKP3hR4phanryHkpPOnbS/qwJKRKad6GA/sfAAsXCYWGSAWHiQWCRMLDxIMK+IYGE5AOGuVvq2bnDrImGIhIlFBhPvC486A39WLgAdLz3GQO1mty4aIRYZ9MIkTKhiCSUnfRKASF83dfdcSSwaJRYNQyRCLBp2P+yRMOVnXkjOu44EoO1/H6Dlyf9wATGCPzufJd+4I/F+x93fI9LRNOrxFr37Y5S8dxUA/XWbqfvPH435t8m66FeJgOpY/xd67POjbpe99PBEQMXCg7T97wNj73PxwYmAGmh4h66X/2fU7fzZeZQzFFB92+2Y+5wIaQ0oY8xpwDXAu4AG4Dpr7a3GmBDwC+BsIALcaK29Op1lE5GdxWJRYoP9RAf6iYX7iQ0O4AsEySipduvDg3RvfM4LFLd+KFwGyD/8ZDKr9gGga8MzdK7/C9GwFzjeNrHwIL6MTBZdeHPie7ffejHh9sZRy1R43FmUnvxpAAYattLw4I1jlj/v4BMSAdX71kt0b3x29OMc7E8+aPpr3xpzn9HBvqR3vqFw8vnx+QPgDyJP0v4AABRfSURBVODz+/Fn5gz7XEZxJf5gCPx+fH4/+Nx2+AME8ksT2wVyCsneZ5m3jd9tn/Tsz8hMbJuzzxEEc4vA29bn8yVex88RgC8jRPGJnxyxjc/t1+cjs2LJ0D73O5JAbpHbzjd8O18wY9gxlZz4yTH/ThMhbQFljFkI3A+cDzwEHAn82RizBTgJMMC+QCHwiDGmxlp7V7rKJzLTxWJRYgN9RPt7iQ70Eh3oI5CdR0ZxFQDhjia6Xnua6ECv226gz70e7Cc62E/Fh75KsMD9UDb+4Wa6Njwz/Ifbk7X4EOZ96kr3neEBGh766ZhlylpgEgEV7miid8sro27nyxg+q48/Mxd/Vg++QAa+QBBfMIQvGMQXyCCYX5LYLphXTO4Bx+ELBCEQHNreeySHRP7hJ5O1+BB8gYBb7/c+4w8QyCtK+u5s5n3mGnz+AL6AFziBoBcsQfzZQ/O4Fh77YQpXfMht4xt/dLh5n/rBuOvjMqv3ofqTqU1eXrD8fSlt58/IpPiEs1P7/solZFYuSWnbrIUHpLTdnkpnDWoJcI+19kHv/TpjzJPA8bjQ+oy1thVoNcZcD6wGFFAyZ8RiscSPXLijmYGGd4j2dRPp6yLa3+MefT3EImEqPjTUzLLj7jUM1L9NtL8XN8bzkIKjz6TsfZ9z+2xvpOV/7h7z+6N9XeAFVCwaTYSTLxjCF8rC7z0HC8sSn/FlZJJ70PH4MzLxeQ9/MBNfMANfRoiQF04AuQccR6hyidtPMOS2CYZcqARDw8qy4Is3pPQ3C1UsovLjl6a0bc5+y1PazucPkDX/XSlvK5MnbQFlrX0KeCr+3hhTAqwE7gaqgQ1Jm28EDk1X2UQmUiwWdcHS0wnRMKHyRW55NELLk/cQ7ekk0ttJtK+baF8nkd5uon3dlJ3+efIPPxmA7k1/p/nPvxnzO8rPvCDx4xgd6CPa3wO4mog/Mxt/KBtfKHt4baOwnMIVH/bWZeEPZbnXGZn4QpkEvYv0AGWnf4Gy07+ALyM07o+wLxCk8qPfSOnvklFUQUZRxa43FPFMSScJY0wh8DDwPPCCt7gnaZMeIGfk50SmSnSgl0h3u/doI2vhQYkeWh0vPELXa08T6elwwdPblbguEapYxIIv/sTtxOen/e9/gEh41O+I9HYlXmeUVJO9z+H4s/LwZ+USyMrFF8rBn5mDP2v4P42qcy7zmrOyxw2TYEEZpaeen9Lx+jOzU9pOZDKlPaCMMfvjrkFtAFYB8X8Jyf8icoAuRCZRLBYj1t9DuKuVSGcLvoxMshYYwDWH1T/4EyLdrUS623e6FlN97hVkLz3MbdvRRN+214et92fm4M8pIFhQnljm8/koPfnT+IIh/Dn5BLzw8WfluQBKCoWcfY4gZ58jUjqOYNL1E5HZJN29+N6DC6dbgMustTGgzxhTh+skUeNtegDDm/xEdlu0v5dwRyMZJfPcRW6g9ZkH6N38EpGuFsKdLcOCJ3vf5VR/4nLANZX11wx1ofUFQwRyCwnkFhHILcQXGrqon3fYyWQvPZxAbiH+7AIC2XmJ7xup8JgPTsahisxK6ezFty/wB+Bya+3PR6y+G7jCGPMykAdcCvwsXWWTmS3c2ULXq38j3N7oHh1NhDuaiPZ1A7DggpsIlc4HYLClhr6tryU+68vIJJhfQiCvZFjPJX92HtWfvopgXrHrchvKGrOXVqh0HpTOm7wDFJmj0lmDugjIB642xiTf43QzsAa4AXgNN/zSWlwtS+awWDRCpLOFwdY6BlvrCbe558HWerKXHELpKecBEOnpGLV3mi8YIlhYRmxg6L6VwqM/SN4h7yGYX0owvwRfKHvU4PH5fGQvOmjyDk5Edimdvfi+AYzX3eci7yFzTKSvm8HmGgabtpN3yEp8AXczYN29/07v2y+P+plA0r0oGUUVFB7zQYKF5QQLygkWlhEsKMOfU7BT+MTvlheR6U9DHUlahbta6d7wDAPNNV4o1RDpbkusz5z3LkLlCwEIFlYQyC0iWFxJRnEVGUVVQ6+9m0/BdUgoPU0TQIvMNgoomXDRgT4GGrcyUL+FgYZ3yCidR+HRZ7p1PZ00/+X2Ydv7giEySueTUToPkmo8ZWd8Cd+ZX05r2UVk+lBAyYTofPkJujetY6DhHcKt9SSPaJC16KBEQGWUVJO/7DRCZQtcKJXNJ1hQ5sb6GkF36YvMbQooSVm4s5X+HW/QX/sm/TvepOyMLyWa2vpr3hgaVdkfJFQ2n1DlEkIVi8ms3i+xD18wg/IPXDAVxReRGUYBJWOK9vfS/o8/eYH0BpHOlmHr+2o2JQIq79ATyVxgyKxcQkbpvERHBxGRPaWAEmKxGOHWOnq3biDa30PRig8Bbpy11r/9J0Td0Dy+zBwyq/cls3pfsua9i6yFByb2kbXAJEZhEBGZCAqoOSgWizLQsJW+ba/Tt3UDfVs3JHrS+TJzKDz6A26qgWAGJSd9kkBeEZnV+7ma0SjXikREJoMCag5wM5EOJObG6XzxUZoe+fWwbfw5BWQtPJDsRQcRi4QTHRSKjjsr7eUVEQEF1Kw12FpH7+aX6HnrJXrfeYXCo86g5L2fAiBzwQEECsrIXnQQWQsPJGvRQWSUzt/lhGsiIumkgJpFere8QvfG5+jZ/BLh1rph6wbbGhKvQxWLWfzVW9NdPBGR3aKAmqFisRgD9VvIKKrAn5ULQNerT9H5z8cB8Gflkr3kMLL3OYKcfY8gWJA0C6pqSiIyAyigZpBYJEzv1tfo2fg83Zv+TqSrlYqPfJ28Q1YCkHfISgIFpeTscwSZ8/bTja4iMqMpoGaA7k3r6LbP0bPpH0T7huZxDOSXEA0nzWe05FCylxw6FUUUEZlwCqhpKNrfO2z+obZn7qd/xxsAZJQtINccS+4BKwhVLlVznYjMWgqoaSLa103Xxufosc/T8/Y/WfC5awlVLAag4KgzCHccQ645hlDZgikuqYhIeiigplAsGqF38z/pfOVJejatIxYe8Nb46NvxRiKg8g89ceoKKSIyRRRQUyQWi7L91q8z2LIjsSxr8SHkHXQ8OfsfQzCvaApLJyIy9RRQaRLpbqdrw9PkH3Eq/oxMfD4/mQsOAGLkHXoSeYe+h4zCiqkupojItKGAmkSxyCDdb/yDrpefpOet9RCNEMgpJO/gEwAoO/1z+DKy1NFBRGQUCqhJEO5qpXP9X+h44c9D05n7/OTsdySB/JLEdv5Q9hSVUERk+lNATYK6e3/IQP3bgOsWnn/EqeQdvFLXlUREdoMCai/FwoN0vf6MG3DVu4ZUsPx99Ly1nsKjP0DW4kPUhCcisgcUUHso3NFMx4uP0vnSX4h0t1N47EcoPeU8wAVUwfL3TXEJRURmNgXUburf8SZtzz9M98bnIBoB3OjgmVVLp7hkIiKziwJqN7Q+dZ+bAh3A5yf3gOMoOPoMshYepGY8EZEJpoDahehAH/5QFgDZ+xxB27MPUXDk6RQe/YFhU1iIiMjEUkCNIhaL0ffOq7Q+/Xt8gRDVn/w3ALLm78/ii3+dmDpdREQmjwIqSSwWo/ftl2l7+vf0bXsdcBP/hbvaEl3EFU4iIumhgMILprfW0/r07+mv2QSAPyuPwhUfovCoMxIz1oqISPoooIBoTwf1919HLDyAP6eAohUfpuDI9+PP1EgPIiJTRQEFBHILKXr3R/FlZFGw/H2JThEiIjJ1FFCe4pX/MtVFEBGRJP6pLoCIiMhoFFAiIjItKaBERGRaUkCJiMi0pIASEZFpSQElIiLTkgJKRESmpdlwH1QAoK6ubqrLISIiuynptzswct1sCKhqgFWrVk11OUREZM9VA28lL5gNAbUOWAnUApEpLouIiOyeAC6c1o1c4YvFYukvjoiIyC6ok4SIiExLCigREZmWFFAiIjItKaBERGRaUkCJiMi0pIASEZFpSQElIiLTkgJKRESmpdkwksReMcYcDtwCHAZsBj5nrd3pjuaZyhjzOeBWoD9p8UXW2junqEgTwhhzDPAHa22F9z4E/AI4GzeiyI3W2qunsIh7ZZTjywQ6gYGkzf7XWvu+qSjfnjDGnAZcA7wLaACus9beOlvO3TjHN+PPHYAx5oPAj4CluOO7drLP35wOKO8P+xDwU+A9wMeBR40xi621HVNauImzHLjBWvvdqS7IRDDG+IDPA9ePWHUlYIB9gULgEWNMjbX2rjQXca+Mc3yHAi3W2qr0l2rvGWMWAvcD5+P+zR0J/NkYswU4iRl+7nZxfM3M4HMHYIypBv4L+Ki19k/GmOXAM8aYdcD/YZLO31xv4jsJyLDW/tRaO2itvRd4DThnaos1oY4EXprqQkygK4EvA/8+Yvn5wA+tta3W2i24H/jVaS7bRBjr+Gb6eVwC3GOtfdBaG/VaKZ4Ejmd2nLsljH18M/3cYa2tBcq9cPIDpUAYVzOctPM3p2tQwEHA6yOWbcT93+qMZ4wJ4JouP22MuRHoAX4D/NhaO1MHYbzFWrvGGHNSfIExpgg32OSGpO1m6nnc6fg8y4EKY8zLQCXwN+Dr1tqadBdwT1hrnwKeir83xpTgBnm+m1lw7nZxfO9nBp+7OGttpzEmB2jHZcePgUYm8fzN9RpUHu5HO1kPkDMFZZkM5cA/gDtx7cZn4/7v/MtTWai9Ya3dMcriPO85+VzOyPM4xvEBdAPPAKfgmlN6gQfTVa6JZIwpBB4Gngde8BbP+HMXN+L4HmIWnTugD8gFjgY+B1zsLZ+U8zfXa1DdQPaIZTlA1xSUZcJZa+uAE5MWvWSM+TnuWtsvp6ZUk6Lbe04+l7PmPAJYa7+R/N4Y8w2g0Riz0Fq7bYqKtduMMfvjfrQ3AKsYOmez4tyNPD5rbRSYFecOwDueAeAfxpi1wFHeqkk5f3O9BrUB9380yQ5geHV1xjLGHGyMuXLE4hDu/4JmDWttK1DH8HM5a84jgDHmB8aYA5MWhbznGXMujTHvwdUq/hs421rbN5vO3WjH5y2fDefuRGPMCyMWZwKTev7meg3qCcBnjLkE103y47hrNjO1+j1SG/BNY8x24DZgGfA14CtTWqrJcTdwhdfOnwdcCvxsaos0oQ4DjjLGnOu9/xnwR2tt4xSWKWXGmH2BPwCXW2t/PmL1jD93uzi+GX3uPC8B873a38+AFbjeph/FBdSknL85XYOy1g4AZ+CCqQW4HDhrhv2HMybvIuyHcT1qOnDdYK+y1v7XlBZscqwBXsX1wlyHO9ZbprREE+vzuP9bfRPYgmtm+fRUFmg3XQTkA1cbY7qSHj9mdpy78Y5vpp87rLXtwAeAj+F+K9cCX7DW/pVJPH+aUVdERKalOV2DEhGR6UsBJSIi05ICSkREpiUFlIiITEsKKBERmZYUUCIiMi3N9Rt1ZQ4wxtyBG3F5LFfiRp5+Asi31qZlmB1vMN9ngPOstZv24PMx4EPW2j9MQFnKgdOstfdMwL4uBXKttSNHMRHZLapByVxwMW7E5WrcFCsAxyQtux74X+919yifnyxfA/65J+HkqQb+MkFluRZ3E+ZE+DlwnjcuncgeUw1KZj3vLvh2AGNMmbe40RtMN9nI95PGGJMF/Ctw8p7uY5Ty7w3fRO3IWttvjLkTuAz4zETtV+YeBZQI4M2/lGji85rPPokLEYObtuRTwLdww9R0AP9qrb3b+3w+cANuSpMY8D/AxeNMn/EJoM1a+6r3+SXA27ihqW4EFgCP4aZGuR74EFALXGitfdT7TKKJzxjzJPBX4HDgfcA23JTjv/G23QJcb639xYjvO9Qr8/nxfVprfcaYDOAqXMBkA88BX7PWWm+7ld7xHoobxue33t8j4h3fA8ALxphvzZahwyT91MQnMrZrgK8DxwKLgBdxwXQ07gf4VmNMfC6qtbggOx03xUkMN+X3WP8TeCbwyCjLrwLOBU7DNUe+jGt+PBJYj5twcizfwTX5LcOF26+MMalMM349cB9usNNqb9kPvDKegxsY1AJ/NcYUetfO/tvb/kDgPOCLJNWWvOBtxM2BJLJHFFAiY7vZWvuEtfYl3I9xF3CZV4u4EVezWGqM2QdXIzrXWrvO+3H+NG4a8PePse+jcINrjvRDbx9P42Zo3WCtvclauxG4GVjo1dZG86S19mavfP+KayE5bFcH6XUK6QX6rbV1xphs4BLgAmvtX621G621X8M1k34aKARKgFpr7RZr7WO4WttjI3a9gaH5gkR2m5r4RMb2ZtLrHmCLtTY+unJ8Lp9MYLH32hozbHqxHFytarRedpVAUwrfmbxN8nd2jvLZRGcLa22HV5aMUbbblX2973jUa0aMywKMtbbFGHM1sNYY833g/wH/Ya1dN2I/zUDFHny/CKCAEhnP4Ij30TG2C3rbLsM17SVrGeMzUUbvmJDqd45mYJRl8e8YWa7x/u3H150GNIxY1wFgrb3M677/Edw0DH8xxnzPWvujpG0DzKBJ+WT6UROfyN57HVdTybXWvmmtfRPXoeE6YKyu1nVAeZrKBy68CpPe7zNifXKAvQmEgYqk43kL+D5wjDFmkTHmV8A71trrrLXvBX6M60SSrIw09oyU2Uc1KJG9ZK21xpiHgbuMMRfhOgf8ENe5YuMYH3sB1+MuXdYBq40x/w835fi/MzyUuoBDjDGLrbXvGGN+CfzMGNMPvIG7JnUWcAXu+M4GMMbcABTgalsjm/gOw81ULbJHVIMSmRjn47qi/zfuh7oQNzJD2xjb/xHX2y9dLsddo3oWuAtXG0puPrwDqAJe93r+fRv4L+B2XE/C5cAZ1tq3rLW9wAdxXcxfwnWOeBV34zEAxpiDcTPMTtSNxDIHaUZdkSlgjMnBTf/9fmvti1NcnAnndaKosNZ+fqrLIjOXalAiU8Ba24O7RnXRVJdlonnd1M/FDZ8ksscUUCJT5yfAYWZE3/RZ4CvAnfFRJ0T2lJr4RERkWlINSkREpiUFlIiITEsKKBERmZYUUCIiMi0poEREZFr6/xFdT1eslqPpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(coffee_results.T, label='coffee')\n",
    "plot(milk_results.T, '--', label='milk')\n",
    "\n",
    "decorate(xlabel='Time (minutes)',\n",
    "         ylabel='Temperature (C)',\n",
    "         loc='center left')\n",
    "\n",
    "savefig('figs/chap16-fig01.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what happens when we mix them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "63.10924357749633"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mix_last = mix(coffee, milk)\n",
    "mix_last.temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.6806722650686297"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mix_last.temp - mix_first.temp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function takes `t_add`, which is the time when the milk is added, and returns the final temperature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_and_mix(t_add, t_total):\n",
    "    \"\"\"Simulates two liquids and them mixes them at t_add.\n",
    "    \n",
    "    t_add: time in minutes\n",
    "    t_total: total time to simulate, min\n",
    "    \n",
    "    returns: final temperature\n",
    "    \"\"\"\n",
    "    coffee = make_system(T_init=90, r=r_coffee, volume=300, t_end=t_add)\n",
    "    coffee_results = run_and_set(coffee)\n",
    "    \n",
    "    milk = make_system(T_init=5, r=r_milk, volume=50, t_end=t_add)\n",
    "    milk_results = run_and_set(milk)\n",
    "    \n",
    "    mixture = mix(coffee, milk)\n",
    "    mixture.t_end = t_total - t_add\n",
    "    results = run_and_set(mixture)\n",
    "\n",
    "    return mixture.temp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can try it out with a few values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "61.4285713124277"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run_and_mix(t_add=0, t_total=30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "62.9028090119359"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run_and_mix(t_add=15, t_total=30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "63.10924357749633"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run_and_mix(t_add=30, t_total=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then sweep a range of values for `t_add`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "sweep = SweepSeries()\n",
    "for t_add in linspace(0, 30, 11):\n",
    "    sweep[t_add] = run_and_mix(t_add, 30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the result looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving figure to file figs/chap16-fig02.pdf\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXxV1bXA8V8SEkIIMwTCIJOwEAUFBJwQBByqOFetItXaqrXWoQ6t9qlF7OtzHqqCoq1aFW3VKgqKCgrOgiiIDAtkhgQIQ4DMw73vj30SLyHDDST33Jus7+eTT+6Z9lmHo3dl77PP3nHBYBBjjDEm2sT7HYAxxhhTGUtQxhhjopIlKGOMMVHJEpQxxpio1MTvAKKJiDQFhgKZQKnP4RhjTGOQAKQDC1S1MHSDJah9DQU+9TsIY4xphEYAn4WusAS1r0yAl19+mU6dOvkdizHGNHhbtmxh/Pjx4H3/hrIEta9SgE6dOtG1a1e/YzHGmMZkv8cq1knCGGNMVLIEZYwxJipZgjLGGBOVLEEZY4yJSpagjDHGRCVLUMYYY6KSdTM3xhgTtuLSYnKL88kryiO3OJ/CkiL6tOtJ0yZJdX4uS1DGGNOIVEwwuUX55BXnkVuUT25xHnnF+eR62/LKf7ttucX5FJcW71fmqYeO5NdDflHnsVqCMsaYGBMIBsgpzCW7YA/ZBXvIKco9qARTGwlx8TRPSqF5YgopSc1ITWrOsd0G19GV7csSlKkTBQUF5OTk0L59e79DMSYmBYNB8orzy5PObu/3Psv53ufCvQSCgQM6T8UEE/q7eVIzUhJDP7vfofskJSQSFxdXx1dfOUtQDdQ///lPJk+eTEJCAhMmTGDNmjU8/PDDB1XmhAkTGDNmDJdffvl+2y699FJ++9vfMnbs2IM6hzENTUFJoUsy+XvYXRiSZCpJQMWBkrDLbZ6UQuvklrRObkmLpqlRmWAOliWoBurVV1/l9ttv5/zzz4/I+Xbu3BmR8xgTLUoCpWTl7mBrThbb83ZVmXQKSgprLsyT3KRpedJp5f3eb7lZS1o1bUFiQmI9Xl10sATVAJ166qls2LCBu+++mwULFtClSxeWL1/O5MmTefzxx1m3bh15eXl89dVXpKenc9NNN5XXfN5//32mTp3Kxo0bCQaDjB07lkmTJpGYWPX/DNdeey0ZGRncdNNNXHfddVx55ZXMmTOHxx57jIyMDA499FDuuusu+vfvD4CIcM899zBlyhSys7O59NJLGTBgAPfddx/Z2dlceOGF/OlPfyrf99Zbb+X555+nsLCQ8847j1tvvZUmTew/XVP/CooL2JKzna25WWzNyXKfc7LKk1I4zWyJ8U1+SjLNWtG6aQtaN9s/CbVKbklyk6YRuKrYYf+XH4T/++RJvsv8ISLnGpR+BLefeG1Y+77//vuMHj2aP//5z4wdO5bHH398n+3vvfcekydP5u9//ztTpkxh4sSJjBkzhoyMDP74xz/yj3/8g6OPPpr169dz0UUXMXv2bH72s59Veb4nn3xyn/MtWbKEW265hcmTJzN06FBmzpzJFVdcwQcffEDLli3LY5g5cyY//vgjF1xwASeeeCLTp09n/fr1XHDBBZxzzjmICACzZ89m+vTpFBQU8Jvf/IZ27dpx1VVXHeC/pDE/CQaD7C3MYUtOFlu8xLPVS0Jbcrezu2BPlcfGEUf7lLZ0TG1Ph5R25Uln359WNEtMjpkmtWhjCaoRGjBgAKNGjQLgrLPO4sknnyQ3N5cOHTowY8YMunXrxt69e9m5cydt2rRh27ZttSr/9ddf56yzzuLYY48F4Oyzz2batGnMmjWLCy+8EIArr7ySlJQUBg4cSEpKChdddBGpqakcfvjhpKWlsWnTpvIEddNNN9GuXTsArrjiCl588UVLUCZsgUCAHfm7yhNQaC1oa8528ksKqjy2SXwTOjZvT8fU9nRM7UDH1PZ0Su1Ax9QOpDVv1yia2fxkCeoghFujiTZlX/ZAeVNZIBAgMTGRN954g9dff53k5GT69+9PYWEhwWCwVuVnZGTw9ddfM3PmzPJ1JSUlZGRklC+3bt26/HNCQgItWrQoX46Pj9/nnN27dy//3LFjR7Zv316reEzDV1RazLbc7WzZG1ILynW1om25OygN7DfVULmUxGblCahTagcvIbnPbZu1Jj7eBtzxiyUoU27mzJm88847vPHGG3Ts2BFwNazaSktL47LLLuPmm28uX7du3bp9uqDXpslj69at5fFkZGTYbMeNXEFJIet2bWLNrvWs3rmeNTs3kLF3K0Gq/kOqTXKrkFpQBzqFfG6R1Nya4KJURBOUiKQDU4CTgAJgqqreKSKdvPWjgCAwA7hOVXdXUkYccA9wFZAEPAfcqqol3vaTgMeA3sBiYIKqrq7nS2sQ9u7dS0JCAklJSRQXF/Pqq6+iqhQX1/xiX2JiInv37gXgnHPO4frrr2fs2LEMHDiQr776iquvvpqpU6dyzDHH1Dquxx57jEcffZTs7Gyee+45Lr300lqXYWJTUWkx67M3lSei1bvWs2lP5n61+vi4eNKat3eJp3mHfZrj0lLbW+eDGBXpGtR0YCHQEUgH5onIcuASYBvQGUgE3sAloesrKeMq4DxgMFAIvAn8GZgkIu2Bt4ArgLeBG4H3RaSvqh7YW22NyLnnnsvXX3/N2LFjSUpKYvDgwYwbN45Vq1bVeOx5553H3Xffzdq1a7npppu44447uOOOO9i0aRNpaWlMmjTpgJITQI8ePTjzzDMpLCxk/PjxjB8//oDKMdGtpLSEDbs3s9pLRGt2rmfj7gxKK/SUi4+Lp3vrLvRq253ebbrTu+0hdGvV2Z4HNUBxtX2+cKBEZDgwE0hX1WJvXS8gH9gBoKpFItIN+DcwS1UnVVLOF8DzqjrVWx4LvKCqXUTkKuByVT0uZP+NwBWq+mEYMfYA1s6ZM4euXbse3AWbOiEivPXWWxx22GF+h2LqUGmglE17Mlm9cwNrdq5n9a71rM/eTEmFF1Xj4uLo2qKTS0beT/dWXUiqh4FJjT82bdrEmDFjAHqq6rrQbZGsQQ0BlgATReRyXBPfZFV9qGwHEXkFuAhYBkyuopz+3vYyK4DOItK2km0ACgwAakxQxpi6FwgEyNi7ldVeIlqzcwNrszdWOiZc5xYdvZrRIfRu250erbuSnJjsQ9QmGkQyQbUFRgDzgF5AP2CWiGSq6jRvn18B1wD/BP4LnFhJOalAXshy2eeUSraVbU+piwswxlQvEAywJSfL1Yp2bmDNrvWs2bWRwkpGU+iY2oHebQ4prx31bNONlMRmPkRtolUkE1QhsEdVJ3rLi0XkWdzzpGkAqloAFIjIH4FVItJWVSuOoZMLhP5XXJZ8crxtFZNRirfNxCBV9TsEU42SQCkrsn5k0ZZlrN65jjW7NpBfvP97RR1S2pYnol5tDqFXm0NIbdrch4hNLIlkgloBpIhIkqoWhZ5fRBYAf1TVj731TYESXMKpaBkgwOfecj8gU1WzRWQZUPEJej/g/rq7DGMat5yiXBZlLuWbjCUsylxKXnH+PtvbNmu9TzNdrzaH0DK5RRWlGVO1SCaoD4Es4CERuRmXZH6Na9I7GbhbRBYDccCDwL9UtbJRFl8EbhGRObgENtFbB65H3/0icqH3+QYgAMytp2syplHI3LuNbzZ/z8KM71mxffU+Y9B1admJIZ0HcliHQ+nV5hDaNGvlY6SmIYlYglLVAhEZCTwOZOI6Sdyvqm+IyCzgIVwtqwR4Dbi97FgRyQGuVtWXgadw3dS/wDXfvQbc5Z1jm4iciXsP6h/AUuDMkBqbMSYMpYFSdPsaFmZ8z8KMJWTs3Vq+LSEuniPShCGdBzCky0A6pXbwMVLTkEWsm3kssG7mpjHLK8pn0Zafmu5yin5qYW+elMKg9CM4uvMAjuzUn+ZJ1u/I1I1o6WZujIkyW3OyWJixhIUZ37Ns26p9XopNb5HGkM4DObrzAKR9bxLiE3yM1DRGlqCMaUQCgQCrdq71nictYdOezPJt8XHx9O/QhyGdBzKkywA6t+joY6TGWIIypsHLLy5g8ZZlLMxYwreZP7C38Ke3LlISm3FUp/4M6TyQQemHW9dvE1UsQRnTAGXl7ihvulu6bdU+Qwh1TO3AkM4DOLrzAPp16EMTa7ozUcoSlDENQCAYYPXO9a7X3eYlrN+9uXxbXFwc0r63l5QG0qVlJ5tewsQES1DGxLC9hTl8uPpT3v9xHrvyf5qdplmTZI7s1J8hnQcwqPMRtGya6mOUxhwYS1DGxKCMvVt5Vz9i7rovKfIGXe2Q0ra8g0P/Dn1s+gkT8yxBGRMjgsEgy7NW8Y7O5tuMH8pnkB2Ufjhn9B3DgI79rOnONCiWoIyJciWBUr7csJAZK2ezdtdGABLjmzCix3DG9R1D11bpPkdoTP2wBGVMlMopymXO6s95b9XH7MzPBqBl01ROPXQkpxx6Iq2SW/ocoTH1yxKUMVFma04W7678mI/WflE+j1KXlp0Y13cMI7oPs9lkTaNhCcqYKBAMBtHta5ixcjYLNi+mbIzMAR37MU7GcGSn/sTHxfscpTGRZQnKGB+VBkr5etMiZupsVu1cB0BCfAIndB/KOBlD99Y2aLFpvCxBGeODvOJ8PlrzBe+t/IisPDdpdGpSc045dASnHjrK5lQyBktQxkRUVu4O3lv5MXPWfE5+iZsaPT01jTNkNCN7HEtTe75kTDlLUMZEwI871jFDZ/PVpu/KZ6Pt36EP42QsgzsfYc+XjKmEJShj6kkgEOCbjO+ZobNZsX014GajPaH7MMb1HU2vtt19jtCY6BZ2ghKReKAfkAaUAluAH1XVpuQ1JkRBcQEfr/2Sd1d+xNbc7YCb1mJs7xH8rM8o2qW08TlCY2JDjQlKRE4EbgDGAi1CNgWBXSIyC5isql/UT4jGxIYdebuYtWous1d/Sm5xPgBpzdtxet/RjO55HMmJyT5HaExsqTJBiUgf4GngEOBN4DxgGbADiAc6AEcCJwKvishq4GpVXVnfQRsTTXIKc3lp8X+Zt+6r8inTpX1vxskYhnY+kvh4e75kzIGorgb1EjBJVWdWsX2j9zNDRP4EnOMdM6xuQzQmes3ftIhnFr7C7oI9xMXFcWy3IYyTMfRp19Pv0IyJedUlqGPCfb7k7femiLxV3X4ikg5MAU4CCoCpqnqniKQBjwFjgDjgPeAGVd1VSRlLgdCnywlAMnC8qn4hIpOA24CikH3GqerccK7FmHDsKdjLP7/9N19sXAjAYR0O5eqhl9K5RUefIzOm4agyQZUlJxFpC2SraqBsm4gMBtar6o7KjqnGdGAh0BFIB+aJyHLgF8BuoCeQCLwIPAlcUklch4cui8irQHHIM7DBwPWq+lQNsRhzQL7cuJB/LHyVPYU5NE1IYvyR53LKoSdaV3Fj6li1/0eJyM3ABmB4hU0PAptE5PfhnkhEhgO9cMmjQFXXAqOAeUAAuFtVc1U1G3gGOCGMMsd7sf02ZPUQYFG4cRkTrt0Fe3j482d45Itn2VOYw+FpfXnwtDs4rc8oS07G1IPqOklcDPwNuBP4vsLmc3BJ4UERyVTVN8I41xBgCTBRRC7HNfFNVtWHvPIqlv9ddYWJSDNcorxaVXO9delAJ+A2ETkW16HjAVV9Loz4jKlUMBjk8w3f8Ny3/2ZvUS7JTZpy6ZHnMbb3CZaYjKlH1T2D+gPwR1V9rOIGVd0D3C8iicAtQDgJqi0wAldj6oV7p2qWl+Cmle0kIrfgEtRxNZR3BZChqm+HrEvzyn8CuMAr420R2VZNZw9jqpSdv5tnFr7Cgs2LATe6+G+HXkqH5u18jsyYhq+6BNUPmFHD8W8AfwzzXIXAHlWd6C0vFpFncd3Xp3nJ7nHgTGC0qq6oobxf4zpclFPVxbhmwzLzRORF7xyWoEzYgsEgn66fz3Pf/YfcojyaNUnml0edz+hex9u06sZESHUJKg9IDeP4gjDPtQJIEZEkVS3rYdcEQERaAO/gXgQepqqbqytIRLoBA4HXKqw/AThaVR8NWZ1UixiNYWdeNlMXTuPbjCUAHNWpP1cNHU/7lLY+R2ZM41JdgvoEuBRYXM0+lwHfhnmuD4Es4CGv84XgakHXAK/iOmyMUNW8MMo6BljmdagIlQ/cKyKrcF3VR+N6Ao4OM0bTiAWDQeau/ZIXFr1OXnE+KYnNuHzQBYzscYzVmozxQXUJ6v+Az0RkN/BIWUcEABFJBW4GrgNOCedEqlogIiNxzXiZuFrN/cAq4HRcE+A2ESk7JFtVu3rny8F1hnjZ29YDyKjkHAtFZIJX7n9wPRAvV9X54cRoGq/teTuZuuBlFm1ZBsDgzgO4asgltE1p7XNkxjRecWVTS1dGRE4HngNaAYp7V6kNrvazDbhWVadHIM6IEJEewNo5c+bQtavNZNoYBINB5qz5nBcXvUF+SQHNk1L41aALGdF9mNWajImATZs2MWbMGICeqroudFu1g8Wq6rsi0gs4CzfuXhtgO/A18IGq2rMdE7Oycnfw9IKX+X7rcgCGdjmS3wy52GazNSZKVPceVIKqlnpNe694P9UqO6YuAzSmrgWCAWav/pSXFr9JQUkhLZKac8WQiziu29FWazImilRXg1ogIvcCr9U0hJGINMENV3QzMKgO4zOmTm3NyeKpBS+xdJsbdP+YroO5YshFtE5u6XNkxpiKqktQZ+M6NDwmItOBWcBSXBNfHD9NtzES+Dkwn/1HhDAmKgSCAd5fNY9p379FYWkRLZum8ushv+DYbkP8Ds0YU4XqBovdCJzjDQx7HfAUbqSG0NrUVlx37jNU9Zv6DNSYA7Vl7zamLHiJ5VmrADjukKO5YtCFtExuUcORxhg/1Tijrqp+C/wKQEQOwY1EHgC21PRCrTF+CgQCvLvqY15dMp2i0mJaJbfkyiEXM6zrUX6HZowJQ40JKpSqbsC9W2RMVMvYs4Up819Ed6wBYET3YVw+6AJaNK1pcBRjTLSoVYIyJtoFAgFmrJzNv3+YQXFpMW2SW3Hl0ZdwdJeBfodmjKklS1Cmwdi0O5Mp8//Fqp3rABjV41h+Oeh8UpOa+xuYMeaAWIIyDcKsVXP516I3KAmU0LZZa64eOp5B6Uf4HZYx5iDUKkF5Y/D1AZYBSaq6t16iMqYW3lr+PtO+fwuA0b2O55dHnk9KUjOfozLGHKywEpSIJAGPAld5q/oC93mz2o5X1d31FJ8x1Xp96bv854d3iCOOq4deyuheNc1zaYyJFeHOV30PbnbaEfw0t9IDuFHFH677sIypXjAY5D8/vOOSU1wc1w6/zJKTMQ1MuAnqQuD3qvol3ou63hQWV+JmwDUmYoLBIK8smc7rS98lPi6e64b/ihN7DPc7LGNMHQv3GVQasKWS9XuAlLoLx5jqBYNBXlz8X2bobBLi4rn+2CtsuCJjGqhwa1CfADeELAe951J3Ap/VeVTGVCIYDPLCd6+55BSfwB+Ou9KSkzENWLg1qBuAWSJyMpAMPI/rzVdKmDPqGnMwAsEA//z233zw4yc0iW/CTcddaS/fGtPAhZugVgGHARcDh3vHvQy8pKp59RSbMYBLTs988wpz1nxGYnwTbjnhanvHyZhGINwE9S1wmao+X4+xGLOfQCDAUwteYu66L0lMSORPJ1zDwE6H+R2WMSYCwk1Q6UBRfQZiTEWlgVImz/8Xn66fT9OEJP404ncc0VH8DssYEyHhJqipwAwReRpYA+SHblTVd+s6MNO4lQZKefzr5/liwzc0bdKU20dcS/+0Pn6HZYyJoHAT1B3e7/sq2RYEEsIpRETSgSnASbgXfqeq6p0ikgY8BozBzdb7HnCDqu6qopwNQDt+mjxxs6qKt+0kr6zewGJggqquDic+Ex1KAqU89uU/+HrTdzRrksyfR/4ead/b77CMMREWVoJS1XC7o9dkOrAQN+lhOjBPRJYDvwB2Az2BROBF4EngkooFiEh7oAvQUlVzK9n2FnAF8DZwI/C+iPRV1UAdXYOpRyWlJTzy5bMs2LyYlMRm/M/I6+jTrqffYRljfBDuWHzVvowbTk8+ERkO9AKOV9ViYK2IjAIKcSNV3F2WcETkGeCJKooaAqyqmJw85wFLVfUNb/kBEbkeVzP7sKYYjb+KSot5+POpfJv5A82TUrhz5PX0atvd77CMMT4Jt4kvh5+a0yoTThPfEGAJMFFELsc18U1W1YeAcyrsew7wXRXlDAbiRWQ+rsb1LXCjqi4H+uNGWg+lwAAsQUW1opIiHvj8aRZvWUaLpObcOeoGerTp5ndYxhgfhZugTqrkuN7AzcCfwiyjLW6w2Xm4mlQ/3Mu/mao6rWwnEbkFl6CqGvmzFJgP3AbsBO4C3hOR/kAqULE2l4cNxxTVCkuKuP+zKSzZuoKWTVO5a9SNHNK6i99hGWN8Fu4zqHmVrJ4jIj/iRjV/K4xiCoE9qjrRW14sIs/imuWmiUgi8Dhu8NnRqrqiiljuD10WkduB3+FqaLnsn4xScDVAE4UKigu477MpLN22klbJLfnLqBvp2ird77CMMVHgYDs/ZOCa1cKxAkjxxvAr0wRARFrgmuCGAsNUdVFVhYjIjSJyQsiqBK+cAlzzXsUXZfqxf7OfiQL5xQX87ZMnWLptJW2SWzHxpD9YcjLGlAu3k8TplaxuhWviWxzmuT4EsoCHRORmXCL5NXAN8CouWY4Io8NFD2CCiIwDsnFd31fhnkWtB+4XkQuBN3FjCAaAuWHGaCIkryifv33yBCt3rKFdszbcddKNpLdI8zssY0wUCfcZ1IxK1hUBC4Dfh1OAqhaIyEhcM14mrsZzPy65nI5rAtwmUl4BylbVrgAikgNcraov4549PYTrRNEc90zrTFUt9Y4/E/ce1D+Apd42GwUjiuQU5fK/8x5n9c71tE9py19OupGOqR38DssYE2XigsHqOuc1LiLSA1g7Z84cunbt6nc4DdLewhz+Ou/vrN21kbTm7fjLSX+gQ/N2fodljPHJpk2bGDNmDEBPVV0Xui2sZ1AiskZE2layvrOIbKuTKE2Dt6dgL5PmPsbaXRvplNqBiaNvsuRkjKlSlU18InIWUNYZoQcwSUQqPh86tJ7iMg3M7oI9TJr7GBt3Z5DeIo2/jPoDbVNa+x2WMSaKVfcMajFuqKA4b3kQ+45oHsR1376sfkIzDcWu/N1Mmvsom/dsoWvLdO4adQOtm7XyOyxjTJSrMkGp6npgNICIPIcbvHVPpAIzDcPOvGzunvsImXu3cUirLtw56npaJbf0OyxjTAwI90XdX4lIExHpwk/DGsUBTYEhqvpKfQVoYtf23J3cPfdRtuZk0aN1V+4YdQMtm6b6HZYxJkaE+x7UOOA53HBFFe0CLEGZfWzL3cHdHz9CVu4OerU5hDtGXk9q0+Z+h2WMiSHhjiTxf8AHwDBgL25svouBLYT5HpRpPLbkZDHxo4fJyt3BoW17cOeoGyw5GWNqLdwE1QeYpKoLcSM2NFfV/+CS0631FZyJPZl7tzHxo4fZnrcTadeLO0ZeT/MkG6vXGFN74SaofNyQQQArgSO9zwuBvnUdlIlNm/dsYeJHD7MzP5vDOhzKn0deR0pSM7/DMsbEqHAT1CfAnSLSGvgGONcbfXwUYD37jEtOHz/CroLdHJ7Wl9tP/D3NEpP9DssYE8PCTVA340YavwKYhhsDbw+u48Tj9ROaiRVFpcU8/MUz7C7Yw8COh3HbiGtJbtLU77CMMTEu3ASVoKqHAVO80caHAecDx6nqvfUWnYkJ/17ythshIjWNW064mqZNkmo+yBhjahDuaObzRGScqn4D4CWpd+svLBMrlm1bxQydQ3xcPL8/5nKrORlj6ky4NaidgA2cZvaRX1zAk/NfIEiQcw87jT7tevodkjGmAQm3BvUxMENEPgJW43r1lVPVP9Z1YCb6vbDodbJyd9CzTTfO7/8zv8MxxjQw4Sao/sCXQDPgiArbbEKpRuibzd/z0ZrPSYxvwu+HX06ThHD/UzLGmPCEOxbfSfUdiIkdewr28vSClwC4eODZdGvV2eeIjDENUdh/9opId+Ba3Iu51wCnAitU9at6is1EoWAwyNSF09hduJf+Hfpwet/RfodkjGmgwp1RdziwFDeCxM9wTX1HAZ+IyJn1F56JNp+un8/8TYto1iSZ3w2/jPi4cPvZGGNM7YT77fIA8FdVPRVv0kJVvRGYBNxTT7GZKLM9byf/+PZVAC4fdAFpNl27MaYehdvENwj4VSXrXwb+HO7JRCQdmIIbDb0AmKqqd4pIGvAYMAY3z9R7uAkSd1VSRgrwMHA2bj6qT4HrVHWDt30ScBv7zv47TlXnhhun2V8gGGDK/H+RX1zA0V2OZFTPY/0OyRjTwIWboHbgRjRfXWH9UGBrLc43HTfAbEcgHfcC8HLgF8BuoCeQCLwIPAlcUkkZ9wKH4noT5uES26vAcd72wcD1qvpULeIyNXh/1TyWbFVaNk3l6qMvIS4uzu+QjDENXLhNfE8AT4vIxbgazlEicj0wGXg6nAK851i9cMmjQFXX4gabnYcbKf1uVc1V1WzgGeCEKopKBiaq6g5VzfdiGy4iZcl2CLAozOsyYdi8Zwsvff8mAFcdPd6mbDfGRES43cwfFJG9uIkLU4DXcZMV/hVXgwnHEGAJMFFELsc18U1W1YeAcyrsew7wXRWxXFXJvj+oaonXhNgJuE1EjsXV/B5Q1efCjNFUUBIo5Ymvnqe4tJiRPY5hWNej/A7JGNNIhN3NXFWfxtWimuMGj63tNBttgRG4GlMvoB8wS0QyVXVa2U4icgsu6RxXaSkhROQXwB+B071VaV75TwAXeGW8LSLbVHVmLeM1wFvLZ7F613rap7TlV4Mu9DscY0wjUpv3oHoDv8E9+ykVkUXAs6q6KcwiCoE9qjrRW14sIs8C5wHTvPmlHgfOBEar6opqYokD7gT+AJyjqp8AqOpiXLNhmXki8qJ3DktQtbR653peX+rGBP7dsF/a5IPGmIgK9z2osbj3oE4DMoAsXC1nmYgMC/NcK4AUEQmdi6GJV34L4ENcp4thqlrlMyQvkb0CXAacoKqzQ7adICI3VjgkCdecaGqhqKSIJ756nmJ5DlwAABuLSURBVEAwwOl9R3NER/E7JGNMIxNuDepB4BFVvT10pYjcj3sGFU6f4w9xie0hEbkZEODXuFEpXsUlyxHeVB7VeRT3wvAxqppVYVs+cK+IrMJ1VR+N6wlowx3U0rQl09m8dwtdWnbikgFn+x2OMaYRCjdB9QUqewDxDG74oxqpaoGIjMQ142XiajX3A6twz5AKgW0i5X+pZ6tqVwARyQGuxjXT/RYoAdaG7AvQRVUXisgEr9z/ABuAy1V1fpjXaYAftirvrvyIhLh4rht+OUk2AaExxgdhT1iIS1B/rbD+ZODzcE+mqmuAMyrZVO1LNaqaGrKYUMO+rwGvhRuT2VdeUT5Pzn8BgPMPP51ebbv7HJExprEKN0F9A/yPVwP6FFeDGYwbzeEVr6kPsLmhYt1z3/2HHXm76N22O+ccdprf4RhjGrFwE9QJwFfe/qFTb3wGdPN+wOaGimnzNy1i3rqvSExIdHM8xVdbWTXGmHpl80EZALIL9vD0Ny8DcOnAc+nSspPPERljGrvavAfVCzgMN0BrqKCqvlmnUZmICgaDPL3gJfYW5jCgo3Bqn5F+h2SMMeElKBG5FTdIaxAorrA5iBv+yMSouWu/ZGHGElISm3HNsF/aHE/GmKgQbg3qVuB/gPtVNVCP8ZgI25a7g+e/c50erxh8Ee1T2vockTHGOOH+qRwH/NeSU8MSCAaY/PUL5JcUMLzrIEZ0D3dQEGOMqX/hJqgncSOE2xubDci7Kz9iWdYqWiW35MohF9scT8aYqBJuE99rwCfAJSKyBTd/UzlV7VXXgZn6tXF3Bq98Px2A3w69lJbJLXyOyBhj9hVugnoJWIcbM6+msfJMlCspLXFzPAVKGN3reIZ0HuB3SMYYs59wE1Q/4EhVXVmfwZjIeH3Zu6zN3kiH5u247Kif+x2OMcZUKtxnUAtwo4+bGLdy+xreXD6LOOK4dthlNEtM9jskY4ypVLg1qJeBf4rIK8BqKrwLpaqT6zowU/cKS4p48usXCAaDnClj6Z/Wx++QjDGmSuEmqNuAHNxstxUFAUtQMeClxf8lM2cb3Vp15qIBZ/kdjjHGVCvcsfh61ncgpn4t3rKM93+cR0J8gpvjKSHR75CMMaZatRmLLx435Xtf4Hnv9wpV3VM/oZm6klOUy5T5LwJwweFn0KNNtxqOMMYY/4XVSUJE0oFFuG7mDwJtcc1+y0Wkf/2FZ+rCP7/9Dzvzs+nTridn9zvF73CMMSYs4fbiewxYBnQA8r11l+LmiHq0HuIydeTLjQv5bP18miYk8fvhl5NgczwZY2JEuAlqNDBJVQvLVqhqHm4A2eH1EZg5eLvyd/PMN68AMOGo80hvkeZzRMYYE77aDBZb2Qsz7YGiugvH1JVgMMhTC14ipyiXIzv15+TeJ/odkjHG1Eq4CeoN4CHvWVQQQEQG4gaRnV5PsZmDMGfNZ3yX+QPNk1K4ZugEGwjWGBNzwu3FdxPwDLDZW14GJALvADeHezIvwU0BTgIKgKmqeqeIpOGec43B1dbeA25Q1V2VlBEH3ANcBSQBzwG3qmqJt/0kr6zewGJggqquDjfGhmBLThYvLHoDgN8M+QVtU1r7HJExxtReWDUoVc1R1YuBQ3Ev614EHKaq56jq7lqcbzqQCXQEjgEuE5FLgGeBEqAn0Adog6udVeYq4DxgsLfvUODPACLSHngLuBtoDbwJvO91kW8UAoEAT379AoUlhRzXbQjHHzLU75CMMeaAhDvl+xrgaFVdA6wJWd8ZWKSqNT59F5HhQC/geFUtBtaKyCigELgQuFtVc719nwGeqKKoy4BHVXWTt+9E4AVgEi5xLVXVN7x9HxCR63E1sw/DudZY947ORrevpk1yK34z5GK/wzHGmANWZYISkbOAE7zFHsAkEak41cahtTjXEGAJMFFELsc18U1W1YeAcyrsew7wXRXl9Mc1MZZZAXQWkbaVbANQYACNIEGtz97Ev394B4Brhk0gtWlznyMyxpgDV10NajFwI+6ZEMAg9u2xF8SNz3dZmOdqC4wA5uFqUv2AWSKSqarTynYSkVtwCeq4KspJZd85qco+p1SyrWx7Spgxxqzi0mKe+Op5SgIlnNx7BEelH+53SMYYc1CqTFCquh73/hMi8hyu08LBDGtUCOxR1Yne8mIReRbXLDdNRBKBx3HPuEar6ooqyskFmoUslyWfHG9bxWSU4m1r0D5c/Snrd2+mY2oHJhx5nt/hGGPMQQu3k8Sv6mDMvRVAiogkhaxrAiAiLXBNcEOBYaq6qJpylrHv3FT9gExVza5kW9n2is1+DUpJoJR3dDYAvzzqfJJtjidjTAMQ9mCxdeBDIAv3PtXNuETya+Aa3Bh/8cAIb4SK6rwI3CIic3A1poneOnC99u4XkQu9zzcAAWBunV5JlPl8/QJ25O2iS8tONn27MabBiFj3a1UtAEbinj9lArOA+4FVwOnAMGCbiOR4P5vKjvWWx3uLTwGvAV94xy4D7vLOsQ3XRHg7sBP4OXCmqjbY0S4CwQBvr/gAgLPkZOLjGk2PemNMAxfJGhReN/UzKtlU7TAHqpoa8jkA/MX7qWzfT3AdOhqF7zKXsnFPJm2btWZE92F+h2OMMXWmum7mYfd8C6NZztST6cvfB+CMvmNokhDRvzeMMaZeVfeNloM37l4YbA4HH+j21azYvprmic0Y2/uEmg8wxpgYUl2COiliUZgDMn25e/Z0yqEjaWY994wxDUx170HNC6cAb7gjE2GbdmfyTcb3JMY34Wd97W8JY0zDE+5YfH2BB3BDCZU158UBTYG0cMsxdeftFW7kplE9j6V1ckufozHGmLoXbp/kKcAh3u/OwGTclBjtgSvrJzRTlR15u/h0w3zi4uI4s9/JfodjjDH1ItwEdQxwtao+jBuj7ytV/R3wJ+CS+grOVG7myo8oDZRyTNfBdErt4Hc4xhhTL2oz5fsW7/MKfnrP6C3cvEwmQnKKcpm9+lMAzu53is/RGGNM/Qk3QS3GDeoKsBQ3IgRAF2p4ydbUrQ9+/ISCkkIGdOxHr7aH+B2OMcbUm3A7N0wEpotIMfAScIeIzMUNxPpu/YRmKioqKeK9lR8DVnsyxjR84Y5m/j7QF3hPVTNwczXNBx7BOklEzNx1X7G7cC8923RjQMd+fodjjDH1Kuzu4aq6IeTzD8Af6yUiU6lAIFA+pcbZ/U4lLs5aVo0xDVu470H1AR7ETdueSIXnTqqaVvehmVBfbfqOrTlZdEztwDFdG81YuMaYRizcGtTTQDpwH3CwExeaWgoGg0xf4QaFPVPGEh9vU2oYYxq+cBPUUGCUqi6sz2BM5ZZsXcHaXRtp1bQFo3oc43c4xhgTEeH+Kb4RsNFIfTLdm5DwZ31PIqlJks/RGGNMZIRbg7oTmCIi9+Bmsd1nhlpVXVbXgRlnzc71LNm6guQmTTnl0BP9DscYYyIm3AT1mvf73yHrgrjOEkFsPqh6M90bFHZs7xGkJjX3ORpjjImccBNUz3qNwlRqS04WX236loT4BMb1HeN3OMYYE1FhJShVXV/fgZj9vbPiQ4LBICN6DKNtSmu/wzHGmIiqMkGJyDagv6puF5Esqpn+3d6DqnvZBXuYu/ZLAM6yKTWMMY1QdTWoW4G93udb6uJkIpKOm1PqJKAAmKqqd4ZsTwbmAveq6ltVlLEU6B6yKgHXw/B4Vf1CRCYBt7FvR45xqjq3Lq4hUt5b+THFgRKO7nIkXVum+x2OMcZEXHUJ6kVVDQCo6gt1dL7pwEKgI+7F33kislxVp4nIQGAqMLy6AlT18NBlEXkVKFbVL7xVg4HrVfWpOoo54vKLC/jgx3kAnG21J2NMI1Xde1DFIrJP052InCgiTQ/kRCIyHOiFSx4FqroWGAV87E0pPwfXS3BD1aXsV+Z4XEL7bcjqIcCiA4kxWsxe/Rm5xfn0a98bad/b73CMMcYX1SWoykYjnYGbA+pADAGWABNFZLOIrAbOVdVMIAPopaqPUM2zrlAi0gw3PuANqprrrUsHOgG3ichWEVkmIr86wHh9UVJawsyVcwA4+7BTfY7GGGP8E/Zo5p6DGUK7LTACmIerSfUDZolIpqpOO4DyrgAyVPXtkHVpXvlPABfgpgV5W0S2qerMg4g9Yj5dP5+d+dl0a5nOoPTDaz7AGGMaqNomqINRCOxR1Yne8mIReRY3U++BJKhf4zpclFPVxbhmwzLzRORF7xxRn6ACwQBvey/mntXvFOLjbFBYY0zjFclvwBVAioiEDiZ3QAlSRLoBA/lphIuy9SeIyI0Vdk/C9RiMegszlrB57xbapbTh+O5D/Q7HGGN8VVOCuFxEcirsf6mIbA/dSVUnh3GuD4Es4CERuRkQXC3omlrEW+YYYJmqZldYnw/cKyKrgPeA0cAl3u+oFgwGmb7cDQo7ru8YmsTb6FHGmMatugS1gf2TxxagYqeDIFBjglLVAhEZCTwOZOJqNfer6hs1HeslyatV9WVvVQ9cx4qK51goIhOA+4H/eNdwuarOr+kcflux/UdW7lhD86QUxvQ63u9wjDHGd1UmKFXtUdcnU9U1wBk17LPfeVU1tcLyA8ADVRz/GhWa/mJB2aCwpx06iuREm9nEGGPsKXwU2JC9mW8zlpCUkMjP+ozyOxxjjIkKlqCiwNvqak8n9TyOlsktfI7GGGOigyUon23P3cnn6xcQHxfPmTLW73CMMSZqWILy2YyVcygNBji222DSUtv7HY4xxkQNS1A+yinMZc6azwE4u98pPkdjjDHRxRKUj2b9OI/CkkKO7NSfHm26+R2OMcZEFUtQPiksKeK9VR8DVnsyxpjKWILyycdrv2BvYQ6923bn8LS+fodjjDFRxxKUD0oDpbyjswFXe4qLO5hB4o0xpmGyBOWDLzd+S1buDtJT0xjW5Si/wzHGmKhkCSrCgsEg01e4QWHP7Hcy8fF2C4wxpjL27Rhhi7csZ332Jlont+TEHsP9DscYY6KWJagIm77ifQBO7zuapIREn6MxxpjoZQkqgn7csY6l21bSLDGZU3qf6Hc4xhgT1SxBRVDZs6eTe59ISlIzn6MxxpjoZgkqQjL2bmX+pkU0iW/CGX2jfoJfY4zxnSWoCHlnxWyCBDmxx3DaNGvldzjGGBP1LEFFwK783cxb9xVxxHFWv5P9DscYY2KCJagIeHflR5QEShja9Ug6t+jodzjGGBMTLEHVs7yifD5Y/QkA5/Q71edojDEmdliCqmcfrv6U/OICDk/ry6HtevgdjjHGxIwmkTyZiKQDU4CTgAJgqqreGbI9GZgL3Kuqb1VTzgagHRD0Vm1WVfG2nQQ8BvQGFgMTVHV13V9NzYpLi3l35UeATalhjDG1Feka1HQgE+gIHANcJiKXAIjIQFxyqnb8HxFpD3QB0lQ11fuRkG1vAXcDrYE3gfdFxJea4ifrvmZXwW66t+rCkZ36+xGCMcbErIjVoERkONALOF5Vi4G1IjIKyBeRvsAc4G9Aeg1FDQFWqWpuJdvOA5aq6hve8gMicj0wBviwDi4jbIFggLfVnfIsm1LDGGNqLZI1iyHAEmCiiGwWkdXAuaqaCWQAvVT1EX5qtqvKYCBeROaLSJaIvC8ih3nb+gPLKuyvwIC6u4zwfLP5ezL3bqNDSluOO2RIpE9vjDExL5LPoNoCI4B5uJpUP2CWiGSq6rRalFMKzAduA3YCdwHviUh/IBXIq7B/HpBykLHXSjAYZPpyNyjsOBlLQnxCJE9vjDENQiQTVCGwR1UnesuLReRZXLNc2AlKVe8PXRaR24Hf4WpoueyfjFKAnAOM+YAsz1rFqp3raJHUnJN6HRfJUxtjTIMRySa+FUCKiCSFrKt1ghSRG0XkhJBVCV45BbjmPalwSD/2b/arV2WDwp7WZxTJTZpG8tTGGNNgRLIG9SGQBTwkIjfjEsmvgWtqWU4PYIKIjAOygfuAVcC3wHrgfhG5ENeD7wYggOsdGBHrszfxXeZSmiYkcVqfUZE6rTHGNDgRq0GpagEwEvf8KROYBdwf0uOuSiKSIyLjvcXbgK+A74BtXnlnqmqpqm4DzgRuxz2f+rm3raiur6cq01e4nnujex1Pi6apkTqtMcY0OBF9UVdV1wBn1LBPj0rWpYZ8LgCu9X4qO/4TYNBBBXqAtuXu4IsN3xAfF884GeNHCMYY02DYUEd1aIbOJhAMcPwhR9OheTu/wzHGmJhmCaqOFJYU8dGazwEb1sgYY+pCRJv4GrouLTrRp11PDmndxe9QjDEm5lmCqiNNmyRx36l/9jsMY4xpMKyJzxhjTFSyBGWMMSYqWYIyxhgTlSxBGWOMiUqWoIwxxkQlS1DGGGOikiUoY4wxUcneg9pXAsCWLVv8jsMYYxqFkO/b/WZ2tQS1r3SA8ePH17SfMcaYupUOrA5dYQlqXwtw09Jn4qaWN8YYU78ScMlpQcUNccFgMPLhGGOMMTWwThLGGGOikiUoY4wxUckSlDHGmKhkCcoYY0xUsgRljDEmKlmCMsYYE5UsQRljjIlKlqCMMcZEJRtJoo6IyJHAU8BAYA1wharu92Z0LBGRK4CngcKQ1deq6gs+hXTARGQYMENV07zlJOAJ4Oe4UUMeVtX/8zHEWqvkmpoCe4GikN2+UNVT/IivNkTkZOBeoA+wDXhAVZ+O5ftUzTXF8n0aB/wN6Im7pvvr8z5ZgqoD3s2ZDjwKnAicD3wgIt1VdY+vwR2cwcBDqnqb34EcKBGJA34NPFhh092AAL2BVsAsEdmsqv+KcIi1Vs01DQB2qmqnyEd14ESkG/AGcBnu/6MhwPsisg4YRQzepxquaQexeZ/SgdeBc1X1PREZDHwuIguAC6iH+2RNfHVjFJCoqo+qarGqvgosBS7yN6yDNgRY5HcQB+lu4BrgrxXWXwb8r6ruUtV1uC/7qyMc24Gq6ppi9X71AKap6puqGvBaHuYCxxO796kHVV9TTN4nVc0EOnjJKR5oB5TgaoP1cp+sBlU3+gPLK6xbgfuLNiaJSAKuuXKCiDwM5AHPAvepaiwN4PiUqt4lIqPKVohIa9zglMtC9oul+7XfNXkGA2ki8j3QEfgEuFFVN0c6wNpQ1U+BT8uWRaQtbtDmF4nR+1TDNZ1GDN4nAFXdKyIpwG5c/rgPyKKe7pPVoOpGKu4LPFQekOJDLHWlA/AN8AKuvfnnuL/ar/EzqNpS1YxKVqd6v0PvWczcryquCSAX+BwYg2tuyQfejFRcdUFEWgFvA18DC73VMXmfylS4punE/n0qAJoDQ4ErgBu89XV+n6wGVTdygWYV1qUAOT7EUidUdQswMmTVIhF5HPd8bbI/UdWZXO936D2L6fsFoKo3hS6LyE1Aloh0U9WNPoUVNhHpi/sCXwaM56f7E7P3qeI1qWoAiOn75F1DEfCNiEwFjvY21fl9shpU3ViG+0soVD/2rfLGFBE5XETurrA6CffXU0xT1V3AFva9ZzF9vwBEZJKIHBayKsn7HfX3TEROxNUw3gJ+rqoFsX6fKrsmb31M3icRGSkiCyusbgrU232yGlTd+BiIE5E/4Lpano97fhNL1faKsoGbRWQT8A9gEHA98Htfo6o7LwJ/8Z4DpAK3AI/5G9JBGwgcLSKXeMuPATNVNcvHmGokIr2BGcD/qOrjFTbH5H2q4Zpi8j7hOnZ08Wp8jwHDcb1Jz8UlqDq/T1aDqgOqWgT8DJeYdgL/A5wTA//BVcl7YHsWrifOHlyX2XtU9XVfA6s7dwE/4HpbLsBd31O+RnTwfo37a/ZHYB2uGWaCnwGF6VqgBfB/IpIT8nMfsXufqrummLxPqrobOB04D/c9NxX4jarOo57uk82oa4wxJipZDcoYY0xUsgRljDEmKlmCMsYYE5UsQRljjIlKlqCMMcZEJUtQxhhjopK9qGsaBRF5HjficlXuxo02/THQQlWjbjgdETkCWAL09EaMrtX2gy3f2ycBN47cL1V15QGcIwicqaozwtj3CWBBLM4/ZuqG1aBMY3EDbsTldNz0KADDQtY9CHzhfc6t5HjjXA8sPpDk5EkHPgxz30nAJBFpd4DnMjHOalCmUfDegt8NICLtvdVZ3qC4oSouG4+IJAO3A6MPtIxK/r2r23ebiMwBrgMmHug5TeyyBGWMx5tfqbyJz2uOuhj3pSy46UcuBW7FDU2zB7hdVV/0jm8BPISbmiQIfATcUNX0GCLSx9v/RNxI0KuAP6vq2972DsDTwClAJvBIheNr2l5tPDUdX4lfANmq+oN3fA9gLW5IrIeBrsBs3JQsDwJneuX+TlU/8I4pb+ITkbnAPOBIL4aNuGnRnw0553+Bf4jI/6pqcQ3xmQbGmviMqd69wI3AMcAhwLe4xDQU9+X5tIiUzS81FZfITsVNVRLETfO93x+C3rTt7+BmIz0GOAr3/Oc5ESkb3fo1oBNwAu5L/08Viqlpe03x1HR8RWcAsypZfw9wCXAyrvn0e1xz6RDgO9xEl1X5E67JbxAuuU0RkdCp0GfjZm4dUkNspgGyBGVM9Z5U1Y9VdRFudOocXC1HcbWGZkBPEemFq2FcoqoLvFrGBNzU36dVUm4z3Bf3daq6QlWX42odbYGO3nQMI4ErVXWRqs4mJIGEsb3aeGo6vgpH4wYDreh/vXN8hptFdpmq/l1VVwBPAt282lxl5qrqk96/5+24Vp2BZRu9KSrW8NOcQ6YRsSY+Y6r3Y8jnPGBdyJT3ZfP3NAW6e59VZJ+pwVJwtZh9eq2pap6ITAEuEZGjgb64KdsBEoAjgEJVDU0I80M+17S9fw3xNKvh+Mp0BLZXsr7iv1HoPqH/RnsrOba8s4Wq7vFiTaywzw4grYbYTANkCcqY6lV87hGoYr8m3r6DcE1poXZW3FlEmuMmsyub7vttXO1sboX94kISYlEl5VS1vaZ4xoRTfgUBIK6S9eH+G1WmsnNWPEcCUFqLMk0DYU18xtSN5bi//Jur6o+q+iOug8ADuNpRRaOAPsAIVf2bqs7E1VDAfUF/j6t1HBVyzOCQzzVtrymemo6vzBagQw371If2WO/KRslqUMbUAVVVEXkb+JeIXAtkAf+L6wCxopJDduCm+r7I6802GHjU29ZUVVeIyLu4HmxX45LJ/RXOV9P2KuNR1ezqjq/CQlyPu4gRkVa45tMFkTyviQ5WgzKm7lyG64r+Fu4LtRVwsqpmV9xRVb/Czbx8L7AMNyPpLbiZVst6rF3sbfsIeAX4e4ViatpeUzw1HV/RTFzHikg6AVd7+i7C5zVRwGbUNcaERURScFOUn6aq30bonK/gegXeE4nzmehiNShjTFhUNQ/3DOvaSJxPRNJxNbYnI3E+E30sQRljauMRYKBU6LteT+4A7lDV/XpBmsbBmviMMcZEJatBGWOMiUqWoIwxxkQlS1DGGGOikiUoY4wxUckSlDHGmKj0/4fu0+Ly6g9IAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(sweep, label='final temp', color='C2')\n",
    "decorate(xlabel='Time added (min)',\n",
    "         ylabel='Final temperature (C)')\n",
    "\n",
    "savefig('figs/chap16-fig02.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use the analytic result to compute temperature as a function of time.  The following function is similar to `run_simulation`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_analysis(system):\n",
    "    \"\"\"Computes temperature using the analytic solution.\n",
    "        \n",
    "    system: System object\n",
    "    \n",
    "    returns: TimeFrame\n",
    "    \"\"\"\n",
    "    T_env, r = system.T_env, system.r\n",
    "    \n",
    "    T_init = system.init.T    \n",
    "    ts = linspace(0, system.t_end)\n",
    "    \n",
    "    T_array = T_env + (T_init - T_env) * exp(-r * ts)\n",
    "    \n",
    "    # to be consistent with run_simulation,\n",
    "    # we put the array into a TimeFrame\n",
    "    return TimeFrame(T_array, index=ts, columns=['T'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how we run it.  From the analysis (see `chap16sympy.ipynb`), we have the computed value of `r_coffee2`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>init</th>\n",
       "      <td>T    90\n",
       "dtype: int64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>r</th>\n",
       "      <td>0.0116102</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>volume</th>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>temp</th>\n",
       "      <td>90</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_0</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_end</th>\n",
       "      <td>30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>dt</th>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>T_env</th>\n",
       "      <td>22</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "init      T    90\n",
       "dtype: int64\n",
       "r                    0.0116102\n",
       "volume                     300\n",
       "temp                        90\n",
       "t_0                          0\n",
       "t_end                       30\n",
       "dt                           1\n",
       "T_env                       22\n",
       "dtype: object"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r_coffee2 = 0.011610223142273859\n",
    "coffee2 = make_system(T_init=90, r=r_coffee2, volume=300, t_end=30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "70.0"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results = run_analysis(coffee2)\n",
    "T_final_analysis = get_last_value(results.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can compare to the results from simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "69.99999985860761"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coffee = make_system(T_init=90, r=r_coffee, volume=300, t_end=30)\n",
    "results = run_simulation(coffee, update_func)\n",
    "T_final_simulation = get_last_value(results.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "They are identical except for a small roundoff error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.4139239112864743e-07"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T_final_analysis - T_final_simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "**Exercise:**  Suppose the coffee shop won't let me take milk in a separate container, but I keep a bottle of milk in the refrigerator at my office.  In that case is it better to add the milk at the coffee shop, or wait until I get to the office?\n",
    "\n",
    "Hint: Think about the simplest way to represent the behavior of a refrigerator in this model.  The change you make to test this variation of the problem should be very small!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "## A refrigerator keeps the milk at a constant temperature,\n",
    "## so it is like a container with r = 0.\n",
    "\n",
    "## With r_milk = 0, it is best to add the milk at the beginning."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
